* local projections DID
* lag, single diff
* two lag, single diff
* clean controls may include future treated. 
* reweighted, or other approach? 

* [] brackets indicate internal census variable name that can't be disclosed. 


clear all
import sas pik year age [annual earnings] [county id] unemp_all unemp_instate min_unemp_instate ///
	plantcap100 plantcap80 plantcap60 plantcap40 plantcap20 ///
	using "/projects/users/########/Snapshot2022/AnalysisData/female_pikloop_lpdid.sas7bdat", case(lower)
rename [county id] ctyfips

keep if age >=18 & age <=65

gen employ = 1 - unemp_all
gen employ_in = 1 - unemp_instate
gen employ_min = 1 - min_unemp_instate

sort pik year
egen long wrkid = group(pik)
sort wrkid year
xtset wrkid year

egen long wrkloc = group(pik ctyfips)
xtset wrkloc year

egen long county = group(ctyfips)

gen asin_bestwages = asinh([annual earnings])

bysort wrkid: egen mxcap = max(plantcap20)
/*
bysort wrkid: egen mxcap40 = max(plantcap40)
bysort wrkid: egen mxcap60 = max(plantcap60)
bysort wrkid: egen mxcap80 = max(plantcap80)
bysort wrkid: egen mxcap100 = max(plantcap100)
*/
sort wrkid year

gen grp20 = 0
replace grp20 = 1 if mxcap>10 & !missing(mxcap)
/*
gen grp40 = 0
replace grp40 = 1 if mxcap40>10 & !missing(mxcap40)
gen grp60 = 0
replace grp60 = 1 if mxcap60>10 & !missing(mxcap60)
gen grp80 = 0
replace grp80 = 1 if mxcap80>10 & !missing(mxcap80)
gen grp100 = 0
replace grp100 = 1 if mxcap100>10 & !missing(mxcap100)
*/

gen trt20 = 0
replace trt20 = 1 if grp20==1 & plantcap20>10 & !missing(plantcap20)
/*
gen trt40 = 0
replace trt40 = 1 if grp40==1 & plantcap40>10 & !missing(plantcap40)
gen trt60 = 0
replace trt60 = 1 if grp60==1 & plantcap60>10 & !missing(plantcap60)
gen trt80 = 0
replace trt80 = 1 if grp80==1 & plantcap80>10 & !missing(plantcap80)
gen trt100 = 0
replace trt100 = 1 if grp100==1 & plantcap100>10 & !missing(plantcap100)
*/

gen stfips = substr(ctyfips,1,2)
egen state = group(stfips)
egen styr = group(stfips year)

replace plantcap20 = plantcap20/1000
replace plantcap40 = plantcap40/1000
replace plantcap60 = plantcap60/1000
replace plantcap80 = plantcap80/1000
replace plantcap100 = plantcap100/1000

gen dcap40 = plantcap40 - plantcap20

gen dcap60 = plantcap60 - plantcap40

gen dcap80 = plantcap80 - plantcap60

gen dcap100 = plantcap100 - plantcap80

/*
* gen ever mover dummy
sort wrkid year
xtset wrkid year
gen mover = 0
replace mover = 1 if county!=L.county & L.county!=.
gen moveout = 0
replace moveout = 1 if county!=F.county & F.county!=.
egen moveinout = rowmax(mover moveout)
bysort wrkid (year): egen movepik = max(mover) 
bysort wrkid (year): egen moveoutpik = max(moveout)
egen evermove = rowmax(movepik moveoutpik)
*/

* do within person-county
sort wrkloc year
xtset wrkloc year

* make sure event_date = . if never treated. 

by wrkid (year), sort: gen byte first = trt20==1
by wrkid (year), sort: gen byte first2 = sum(first)==1
gen trtyr = 0 
replace trtyr = year if first2==1
bysort wrkid (year): egen event_date = max(trtyr) 
bysort wrkid (year): replace event_date = . if grp20==0

gen event_time = year - event_date 

global post_window 6
global pre_window 6
/*
global preobs = $pre_window +1
gen t = _n - $preobs
replace t = . if t> $post_window
*/

xtset wrkloc year


/*
Do 4 versions: lag0, with spatial lags (4)
- for discrete and continuous treatment,
- for employment and asin_bestwages if unemp_instate<0.75
*/


* LP-DID Event Study

* forward changes in outcomes, two years prior as base year: 
* gen forward changes - employ
forval j = 0/$post_window {
	gen D`j'y = F`j'.employ_in - L2.employ_in
}
gen Dn1y = L.employ_in - L2.employ_in
forval j = 3/$pre_window {
	gen Dm`j'y = L`j'.employ_in - L2.employ_in
}
* gen forward changes - log earnings
forval j = 0/$post_window {
	gen D`j'ya = F`j'.asin_bestwages - L2.asin_bestwages 
}
gen Dn1ya = L.asin_bestwages  - L2.asin_bestwages 
forval j = 3/$pre_window {
	gen Dm`j'ya = L`j'.asin_bestwages  - L2.asin_bestwages 
}

capture postutil clear
tempfile lpdid_estimates2022
postfile handle str32 vble double b_emp20bscpl0 vnc_emp20bscpl0 ///
	b_emp20cbscpl0 vnc_emp20cbscpl0 ///
	b_wg20bscpl0 vnc_wg20bscpl0 ///
	b_wg20cbscpl0 vnc_wg20cbscpl0 ///
	using `lpdid_estimates2022', replace


* 1. dummy treat, single diff treat, no lag, not treated as of t+h, not reweighted, surrounding capacity controls.
* employ_in
matrix blpdid = J(1,13,.)
matrix blpdidvar = J(1,13,.)
matrix colnames blpdid = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
matrix colnames blpdidvar = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
forval j = 0/$post_window {
	reghdfe D`j'y D.trt20 D.dcap40 D.dcap60 D.dcap80 D.dcap100	///
		if D.trt20==1 | F`j'.trt20==0,  ///
		absorb(styr) vce(cluster wrkloc)
	matrix blpdid[1,`j'+$post_window + 1] = _b[D.trt20]
	matrix blpdidvar[1,`j'+$post_window +1] = e(V)[1,1]
	
	if `j'>2 & `j'<=$pre_window {
		reghdfe Dm`j'y D.trt20 D.dcap40 D.dcap60 D.dcap80 D.dcap100 ///
			if D.trt20==1 | trt20==0, ///
			absorb(styr) vce(cluster wrkloc)
		local p = $post_window +1 -`j'
		matrix blpdid[1,`p'] = _b[D.trt20]
		matrix blpdidvar[1,`p'] = e(V)[1,1]
	}
} 
reghdfe Dn1y 			///
	D.trt20 D.dcap40 D.dcap60 D.dcap80 D.dcap100 /// 
	if D.trt20 == 1 | trt20 ==0, 	///  /* clean controls */
	absorb(styr) vce(cluster wrkloc)	
matrix blpdid[1,$post_window ] = _b[D.trt20]
matrix blpdidvar[1,$post_window ] = e(V)[1,1]

matrix blpdid[1,$post_window -1] = 0
matrix blpdidvar[1,$post_window -1] = 0	

matrix blp_emp20bscpl0 = (blpdid \ blpdidvar)
matrix rownames blp_emp20bscpl0  = b_emp20bscpl0 var_emp20bscpl0
matrix drop blpdid 
matrix drop blpdidvar

matrix B = (blp_emp20bscpl0)

* asin_bestwages
matrix blpdid = J(1,13,.)
matrix blpdidvar = J(1,13,.)
matrix colnames blpdid = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
matrix colnames blpdidvar = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
forval j = 0/$post_window {
	reghdfe D`j'ya D.trt20 	D.dcap40 D.dcap60 D.dcap80 D.dcap100	///
		if D.trt20==1 | F`j'.trt20==0,  ///
		absorb(styr) vce(cluster wrkloc)
	matrix blpdid[1,`j'+$post_window + 1] = _b[D.trt20]
	matrix blpdidvar[1,`j'+$post_window +1] = e(V)[1,1]
	
	if `j'>2 & `j'<=$pre_window {
		reghdfe Dm`j'ya D.trt20 D.dcap40 D.dcap60 D.dcap80 D.dcap100	 ///
			if D.trt20==1 | trt20==0, ///
			absorb(styr) vce(cluster wrkloc)
		local p = $post_window +1 -`j'
		matrix blpdid[1,`p'] = _b[D.trt20]
		matrix blpdidvar[1,`p'] = e(V)[1,1]
	}
} 
reghdfe Dn1ya 			///
	D.trt20  D.dcap40 D.dcap60 D.dcap80 D.dcap100 /// /* lags */
	if D.trt20 == 1 | trt20 ==0, 	///  /* clean controls */
	absorb(styr) vce(cluster wrkloc)	
matrix blpdid[1,$post_window ] = _b[D.trt20]
matrix blpdidvar[1,$post_window ] = e(V)[1,1]

matrix blpdid[1,$post_window -1] = 0
matrix blpdidvar[1,$post_window -1] = 0	

matrix blp_wg20bscpl0 = (blpdid \ blpdidvar)
matrix rownames blp_wg20bscpl0  = b_wg20bscpl0 var_wg20bscpl0
matrix drop blpdid 
matrix drop blpdidvar

matrix B = (B \ blp_wg20bscpl0)


************************
* Continuous treatment *
************************

* employ_in
matrix blpdid = J(1,13,.)
matrix blpdidvar = J(1,13,.)
matrix colnames blpdid = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
matrix colnames blpdidvar = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
forval j = 0/$post_window {
	reghdfe D`j'y D.plantcap20 D.dcap40 D.dcap60 D.dcap80 D.dcap100	///
		if D.plantcap20>=0.01 | F`j'.plantcap20==0,  ///
		absorb(styr) vce(cluster wrkloc)
	matrix blpdid[1,`j'+$post_window + 1] = _b[D.plantcap20]
	matrix blpdidvar[1,`j'+$post_window +1] = e(V)[1,1]
	
	if `j'>2 & `j'<=$pre_window {
		reghdfe Dm`j'y D.plantcap20 D.dcap40 D.dcap60 D.dcap80 D.dcap100 ///
			if D.plantcap20>=0.01 | plantcap20==0, ///
			absorb(styr) vce(cluster wrkloc)
		local p = $post_window +1 -`j'
		matrix blpdid[1,`p'] = _b[D.plantcap20]
		matrix blpdidvar[1,`p'] = e(V)[1,1]
	}
} 
reghdfe Dn1y 			///
	D.plantcap20 D.dcap40 D.dcap60 D.dcap80 D.dcap100 /// 
	if D.plantcap20 >=0.01 | plantcap20 ==0, 	///  /* clean controls */
	absorb(styr) vce(cluster wrkloc)	
matrix blpdid[1,$post_window ] = _b[D.plantcap20]
matrix blpdidvar[1,$post_window ] = e(V)[1,1]

matrix blpdid[1,$post_window -1] = 0
matrix blpdidvar[1,$post_window -1] = 0	

matrix blp_emp20cbscpl0 = (blpdid \ blpdidvar)
matrix rownames blp_emp20cbscpl0  = b_emp20cbscpl0 var_emp20cbscpl0
matrix drop blpdid 
matrix drop blpdidvar

matrix B = (B \ blp_emp20cbscpl0)

* asin_bestwages
matrix blpdid = J(1,13,.)
matrix blpdidvar = J(1,13,.)
matrix colnames blpdid = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
matrix colnames blpdidvar = pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6
forval j = 0/$post_window {
	reghdfe D`j'ya D.plantcap20 	D.dcap40 D.dcap60 D.dcap80 D.dcap100	///
		if D.plantcap20>=0.01 | F`j'.plantcap20==0,  ///
		absorb(styr) vce(cluster wrkloc)
	matrix blpdid[1,`j'+$post_window + 1] = _b[D.plantcap20]
	matrix blpdidvar[1,`j'+$post_window +1] = e(V)[1,1]
	
	if `j'>2 & `j'<=$pre_window {
		reghdfe Dm`j'ya D.plantcap20 D.dcap40 D.dcap60 D.dcap80 D.dcap100	 ///
			if D.plantcap20>=0.01 | plantcap20==0, ///
			absorb(styr) vce(cluster wrkloc)
		local p = $post_window +1 -`j'
		matrix blpdid[1,`p'] = _b[D.plantcap20]
		matrix blpdidvar[1,`p'] = e(V)[1,1]
	}
} 
reghdfe Dn1ya 			///
	D.plantcap20  D.dcap40 D.dcap60 D.dcap80 D.dcap100 /// /* lags */
	if D.plantcap20 >=0.01 | plantcap20 ==0, 	///  /* clean controls */
	absorb(styr) vce(cluster wrkloc)	
matrix blpdid[1,$post_window ] = _b[D.plantcap20]
matrix blpdidvar[1,$post_window ] = e(V)[1,1]

matrix blpdid[1,$post_window -1] = 0
matrix blpdidvar[1,$post_window -1] = 0	

matrix blp_wg20cbscpl0 = (blpdid \ blpdidvar)
matrix rownames blp_wg20cbscpl0  = b_wg20cbscpl0 var_wg20cbscpl0
matrix drop blpdid 
matrix drop blpdidvar

matrix B = (B \ blp_wg20cbscpl0)



foreach v in pre6 pre5 pre4 pre3 pre2 pre1 tau0 tau1 tau2 tau3 tau4 tau5 tau6 {
	post handle ("`v'") (B["b_emp20bscpl0", "`v'"]) (B["var_emp20bscpl0", "`v'"]) ///
	(B["b_emp20cbscpl0", "`v'"]) (B["var_emp20cbscpl0", "`v'"]) ///
	(B["b_wg20bscpl0", "`v'"]) (B["var_wg20bscpl0", "`v'"]) ///
	(B["b_wg20cbscpl0", "`v'"]) (B["var_wg20cbscpl0", "`v'"]) 
}


postclose handle
preserve
use "/projects/users/########/Snapshot2022/Results/lpdid2_strateq2_female2.dta", clear
append using `lpdid_estimates2022'
save "/projects/users/########/Snapshot2022/Results/lpdid2_strateq2_female2.dta", replace
restore

